/* ***
  **********************************************
  Algoritmo de Monte Carlo - Paralelo
  **********************************************
*** */ 

/***************
BIBLIOTECAS
***************/

#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#include <time.h>
#include <math.h>
#define N_INTERACOES 1000000000
#define NUM_THREADS 2000


/***************
FUNCOES
***************/

/* Box Muller State - Substituindo função RANDoMIZE por eficiência, e questões da rand */

struct BoxMullerState
{
        double x1, x2, w, y1, y2;
        int useLast;
        struct drand48_data random;
};

void initBoxMullerState(struct BoxMullerState* state)
{
        state->random.__init = 0;
        state->useLast = 0;
        
        struct timeval now;
	gettimeofday(&now, NULL);
	state->random.__x[0] = now.tv_usec;
}

double boxMullerRandom(struct BoxMullerState* state)
{
        double randomNumber;
		
	    if (state->useLast)
	    {
	            state->y1 = state->y2;
	            state->useLast = 0;
	    }
	    else
	    {
	            do
	            {
	                    drand48_r(&state->random, &state->x1);
	                    state->x1 = 2.0 * state->x1 - 1.0;
	                    drand48_r(&state->random, &state->x2);
	                    state->x2 = 2.0 * state->x2 - 1.0;
	                    state->w = state->x1 * state->x1 + state->x2 * state->x2;
	            }
	            while (state->w >= 1.0);

	            state->w = sqrt((-2.0 * log(state->w)) / state->w);
	            state->y1 = state->x1 * state->w;
	            state->y2 = state->x2 * state->w;
	            state->useLast = 1;
	    }

        return state->y1;
}

/* Função do cálculo do Pi */

void* threadCalcPi (void *id){

	struct BoxMullerState state;	
	initBoxMullerState(&state);

//	int IDThread;
	int i;
	float iteracoes_totais;
	double val1, val2, val1int, val2int;

	//Recebe o ID da Thread
	//IDThread = (int)longthreadid;
	//printf("Entrei: %d\n", IDThread);

	//Alocando um contador para nossa amada Thread
	float *contador = (float *)malloc(sizeof(float));
	*contador = 0;
	
	//Divide o númer de iterações pelas threads
	iteracoes_totais = N_INTERACOES/NUM_THREADS;
	
	unsigned int x = time(NULL);

	for(i = 0; i < iteracoes_totais; i++){
		val1 = boxMullerRandom(&state);
		val2 = boxMullerRandom(&state);
		//printf("(%.4lf, %.4lf)\n",val1,val2);

        if(sqrt(val1*val1+val2*val2) <= 1)
			*contador+=1;
	}

	pthread_exit((void *)contador);	
}


/***************
PRINCIPAL
***************/

int main (void){
	pthread_t threads[NUM_THREADS];
	void *retorno;
	float result, result2=0;
	int i;

	for(i = 0; i<NUM_THREADS; i++){
		pthread_create(&threads[i], NULL, threadCalcPi, NULL);
	}

	for(i = 0; i<NUM_THREADS; i++){
		pthread_join(threads[i], &retorno);
		result = *(float*)retorno;
		result2 += result;
	}

	result = 4*(result2/N_INTERACOES)*2;
	
    printf("%.6f", result);

	return 0;
}

